/*
 * @(#)Marquardt.java        1.0 2000/05/09
 *
 * This file is part of Xfuzzy 3.0, a design environment for fuzzy logic
 * based systems.
 *
 * (c) 2000 IMSE-CNM. The authors may be contacted by the email address:
 *                    xfuzzy-team@imse.cnm.es
 *
 * Xfuzzy is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by
 * the Free Software Foundation.
 *
 * Xfuzzy is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 */

package xfuzzy.xfsl.model.algorithm;

import xfuzzy.xfsl.model.*;
import xfuzzy.lang.*;
import xfuzzy.xfds.XfdsDataSet;

/**
 * Algoritmo de Levenberg-Marquardt
 * 
 * @author Francisco Jos� Moreno Velo
 * 
 * @see "Fletcher, R., Practical Methods of Optimization, John Wiley & Sons, 
 * Ltd., 1986."
 *
 */
public class Marquardt extends XfslGradientBasedAlgorithm {

	//----------------------------------------------------------------------------//
	//                            MIEMBROS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Factor de adici�n al Hessiano para garantizar que sea definido positivo
	 */
	private double mu;
	
	/**
	 * Valor inicial del factor de adici�n al Hessiano
	 */
	private double initmu;
	
	/**
	 * Factor de incremento
	 */
	private double increase;
	
	/**
	 * Factor de decremento
	 */
	private double decrease;

	//----------------------------------------------------------------------------//
	//                                CONSTRUCTOR                                 //
	//----------------------------------------------------------------------------//

	/**
	 * Constructor
	 */
	public Marquardt() {
		this.mu = 1;
		this.initmu = 1;
		this.increase = 10;
		this.decrease = 0.1;
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS P�BLICOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Devuelve el c�digo de identificaci�n del algoritmo
	 */
	public int getCode() {
		return MARQUARDT;
	}

	/**
	 * Actualiza los par�metros de configuraci�n del algoritmo
	 */
	public void setParameters(double[] param) throws XflException {
		if(param.length != 3) throw new XflException(26);
		initmu = test(param[0], POSITIVE);
		increase = test(param[1], INCREASE);
		decrease = test(param[2], DECREASE);
	}

	/**
	 * Obtiene los par�metros de configuraci�n del algoritmo
	 */
	public XfslAlgorithmParam[] getParams() {
		XfslAlgorithmParam[] pp = new XfslAlgorithmParam[3];
		pp[0] = new XfslAlgorithmParam(initmu, 1, POSITIVE, "Initial Hessian Addition");
		pp[1] = new XfslAlgorithmParam(increase, 10,  INCREASE, "Increase Factor");
		pp[2] = new XfslAlgorithmParam(decrease, 0.1, DECREASE, "Decrease Factor");
		return pp;
	}

	/**
	 * Ejecuta una iteraci�n del algoritmo
	 */
	public XfslEvaluation iteration(Specification spec, XfdsDataSet pattern,
			XfslErrorFunction ef) throws XflException {
		XfslEvaluation eval = null;
		Parameter[] param = spec.getAdjustable();
		if(init) { mu = initmu; init = false; } 

		//+++++++++++++++++++++++++
		// crear jacobiano y error
		//+++++++++++++++++++++++++
		int N = pattern.input.length;
		int M = pattern.output[0].length;
		int P = param.length;
		double[][] jacobian = new double[M*N][P];
		double[] error = new double[M*N];

		//++++++++++++++++++++++++++++
		// calcular jacobiano y error
		//++++++++++++++++++++++++++++
		for(int i=0,n=0; i<N; i++) {
			double[] output = computeOutputGradient(spec,pattern.input[i]);
			for(int p=0; p<P; p++)
				for(int o=0; o<M; o++)
					jacobian[n+o][p] = param[p].oderiv[o]/pattern.range[o];
			for(int o=0; o<M; o++)
				error[n+o] = (output[o] - pattern.output[n][o])/pattern.range[o];
			n+=M;
		}  

		//+++++++++++++++++++++++++++++++++
		// calcular el error total inicial
		//+++++++++++++++++++++++++++++++++
		double preverror = 0;
		for(int i=0; i<error.length; i++) preverror += error[i]*error[i];
		preverror /= error.length;

		//+++++++++++++++++++++++++++++++++++++++++
		// bucle de actualizacion con backtracking
		//+++++++++++++++++++++++++++++++++++++++++
		boolean backtracking = true;
		while(backtracking) {
			double[] desp = compute_desp(jacobian,error);
			for(int i=0; i<param.length; i++) param[i].setDesp(desp[i]);
			spec.update();
			eval = evaluate(spec,pattern,preverror);
			if(eval.error > preverror) {
				mu *= increase;
				for(int i=0; i<param.length; i++) param[i].value-=param[i].getPrevDesp();
			}
			else backtracking = false;
		}
		mu *= decrease;

		return eval;
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Evalua el estado del sistema
	 */
	private XfslEvaluation evaluate(Specification spec, XfdsDataSet pattern,
			double le) {
		SystemModule system = spec.getSystemModule();
		double mxae=0, error=0;
		for(int p=0; p<pattern.input.length; p++) {
			double[] output = system.crispInference(pattern.input[p]);
			for(int i=0; i<output.length; i++) {
				double dev = (output[i]-pattern.output[p][i])/pattern.range[i];
				if(dev<0) dev = -dev;
				error += dev*dev/output.length;
				if(dev>mxae) mxae = dev;
			}
		}
		return new XfslEvaluation(error,error,mxae,le,pattern.input.length);
	}

	/**
	 * Calcula el desplazamiento de los p�rametros
	 */
	private double[] compute_desp(double[][] jacobian, double[] error) {
		int N = jacobian[0].length;
		double[][] hessian = new double[N][N];
		for(int i=0; i<hessian.length; i++) 
			for(int j=0; j<hessian.length; j++)
				for(int k=0; k<jacobian.length; k++)
					hessian[i][j] += jacobian[k][i]*jacobian[k][j];
		for(int i=0; i<hessian.length; i++) hessian[i][i] += mu;
		double[][] inverse = inverse(hessian);
		double[] gradient = new double[hessian.length];
		for(int i=0; i<hessian.length; i++) 
			for(int j=0; j<jacobian.length; j++)
				gradient[i] += jacobian[j][i]*error[j];
		double[] desp = new double[hessian.length];
		for(int i=0; i<hessian.length; i++)
			for(int j=0; j<hessian.length; j++)
				desp[i] -= inverse[i][j]*gradient[j];
		return desp;
	}

	/**
	 * Calcula la inversa de una matriz
	 */
	private double[][] inverse(double[][] m) {
		int N = m.length;
		double[][] y = new double[N][N];
		double[] col = new double[N];
		int[] indx = ludcmp(m);

		for(int j=0; j<N; j++) {
			for(int i=0; i<N; i++) col[i] = 0.0;
			col[j] = 1.0;
			lubksb(m,indx,col);
			for(int i=0; i<N; i++) y[i][j] = col[i];
		}
		return y;
	}

	/**
	 * Funci�n auxiliar para el c�lculo de la inversa
	 */
	private static void lubksb(double[][] a, int[] indx, double[] b) {
		int ii=-1;
		int N = a.length;
		for(int i=0; i<N; i++) {
			int ip=indx[i];
			double sum = b[ip];
			b[ip] = b[i];
			if(ii>=0) for(int j=ii; j<i;j++) sum -= a[i][j]*b[j];
			else if(sum != 0.0) ii=i;
			b[i] = sum;
		}
		for(int i=N-1; i>=0; i--) {
			double sum = b[i];
			for(int j=i+1; j<N; j++) sum -= a[i][j]*b[j];
			b[i] = sum/a[i][i];
		}
	}

	/**
	 * Funci�n auxiliar para el c�lculo de la inversa
	 */
	private static int[] ludcmp(double[][] a) {
		int N = a.length;
		int[] indx = new int[N];
		double[] vv = new double[N];

		for(int i=0; i<N; i++) {
			double big = 0;
			for(int j=0; j<N; j++) if(Math.abs(a[i][j])>big) big = Math.abs(a[i][j]);
			vv[i] = 1/big;
		}
		for(int j=0; j<N; j++) {
			for(int i=0; i<j; i++) {
				double sum = a[i][j];
				for(int k=0; k<i; k++) sum -= a[i][k]*a[k][j];
				a[i][j] = sum;
			}
			double big = 0;
			int imax = 0;
			for(int i=j; i<N; i++) {
				double sum=a[i][j];
				for(int k=0; k<j; k++) sum -= a[i][k]*a[k][j];
				a[i][j] = sum;
				double dum = vv[i]*(sum>0?sum:-sum);
				if(dum>=big) { big = dum; imax = i; }
			}
			if(j!=imax) {
				for(int k=0; k<N; k++) {
					double dum=a[imax][k];
					a[imax][k] = a[j][k];
					a[j][k] = dum;
				}
				vv[imax] = vv[j];
			}
			indx[j] = imax;
			if(a[j][j] == 0.0) a[j][j] = 1.0e-20;
			for(int i=j+1; i<N; i++) a[i][j] /= a[j][j];
		}
		return indx;
	}
}
